/*
 * Critical kernel based parallel thinning algorithm from 'A New 3D Parallel Thinning Scheme Based on Critical Kernels' by
 * Bertrand and Couprie. Based on the reference implementation curve-skeleton algorithm by Michel Couprie.
 * Author: Rune Havnung Bakken, Lars Moland Eliassen. 2010-2012.
 */

#include "thinning_kernel.cuh"

#define S_OBJECT       1
#define S_SIMPLE       2
#define S_2M_CRUCIAL   4
#define S_1M_CRUCIAL   8
#define S_0M_CRUCIAL  16
#define S_CURVE       32
#define S_SURF        64
#define S_SELECTED   128
#define S_ANY        255

#define IS_OBJECT(f)     ((f)&S_OBJECT)
#define IS_SIMPLE(f)     ((f)&S_SIMPLE)
#define IS_2M_CRUCIAL(f) ((f)&S_2M_CRUCIAL)
#define IS_1M_CRUCIAL(f) ((f)&S_1M_CRUCIAL)
#define IS_0M_CRUCIAL(f) ((f)&S_0M_CRUCIAL)
#define IS_CURVE(f)      ((f)&S_CURVE)
#define IS_SURF(f)       ((f)&S_SURF)
#define IS_SELECTED(f)   ((f)&S_SELECTED)

#define HAS_LABEL(label,voxel) ((voxel)&label)

#define SET_OBJECT(f)     (f|=S_OBJECT)
#define SET_SIMPLE(f)     (f|=S_SIMPLE)
#define SET_2M_CRUCIAL(f) (f|=S_2M_CRUCIAL)
#define SET_1M_CRUCIAL(f) (f|=S_1M_CRUCIAL)
#define SET_0M_CRUCIAL(f) (f|=S_0M_CRUCIAL)
#define SET_CURVE(f)      (f|=S_CURVE)
#define SET_SURF(f)       (f|=S_SURF)
#define SET_SELECTED(f)   (f|=S_SELECTED)

#define UNSET_OBJECT(f)     (f&=~S_OBJECT)
#define UNSET_SIMPLE(f)     (f&=~S_SIMPLE)
#define UNSET_2M_CRUCIAL(f) (f&=~S_2M_CRUCIAL)
#define UNSET_1M_CRUCIAL(f) (f&=~S_1M_CRUCIAL)
#define UNSET_0M_CRUCIAL(f) (f&=~S_0M_CRUCIAL)
#define UNSET_SELECTED(f)   (f&=~S_SELECTED)

#define SINGLE_INDEX(x,y,z) ((z)*9+(y)*3+(x))

#define SET_BIT(integer,index) ((integer)|=1<<(index))
#define UNSET_BIT(integer,index) ((integer)&=~(1<<(index)))
#define GET_BIT(integer,index) (((integer)>>(index))&1)

/* Table of topological numbers, ref. Malandain et. al */
__constant__ int topoTab[256][2]={
  1,0,  1,1,  1,1,  1,1,  1,1,  1,1,  1,1,  1,1,  /*  0 -  7 */
  1,1,  2,2,  2,2,  2,2,  1,1,  1,1,  1,1,  1,1,  /*  8 - 1f */
  1,1,  2,2,  2,2,  2,2,  1,1,  1,1,  1,1,  1,1,  /* 10 - 17 */
  1,1,  2,2,  2,2,  2,2,  1,1,  1,1,  1,1,  1,1,  /* 18 - 1f */
  1,1,  2,2,  2,2,  2,2,  2,2,  2,2,  2,2,  2,2,  /* 20 - 27 */
  2,2,  3,3,  3,3,  3,3,  2,2,  2,2,  2,2,  2,2,  /* 28 - 2f */
  1,1,  2,2,  2,2,  2,2,  1,1,  1,1,  1,1,  1,1,  /* 30 - 37 */
  1,1,  2,2,  2,2,  2,2,  1,1,  1,1,  1,1,  1,1,  /* 38 - 3f */
  1,1,  1,1,  2,2,  1,1,  2,2,  1,1,  2,2,  1,1,  /* 40 - 47 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* 48 - 4f */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* 50 - 57 */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* 58 - 5f */
  1,1,  1,1,  2,2,  1,1,  2,2,  1,1,  2,2,  1,1,  /* 60 - 67 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* 68 - 6f */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* 70 - 77 */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* 78 - 7f */
  1,1,  1,1,  2,2,  1,1,  2,2,  1,1,  2,2,  1,1,  /* 80 - 87 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* 88 - 8f */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* 90 - 97 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* 98 - 9f */
  2,2,  2,2,  3,3,  2,2,  3,3,  2,2,  3,3,  2,2,  /* a0 - a7 */
  3,3,  3,3,  4,4,  3,3,  3,3,  2,2,  3,3,  2,2,  /* a8 - af */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* b0 - b7 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* b8 - bf */
  1,1,  1,1,  2,2,  1,1,  2,2,  1,1,  2,2,  1,1,  /* c0 - c7 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* c8 - cf */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* d0 - d7 */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* d8 - df */
  1,1,  1,1,  2,2,  1,1,  2,2,  1,1,  2,2,  1,1,  /* e0 - e7 */
  2,2,  2,2,  3,3,  2,2,  2,2,  1,1,  2,2,  1,1,  /* e8 - ef */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1,  /* f0 - f7 */
  1,1,  1,1,  2,2,  1,1,  1,1,  0,1,  1,1,  0,1   /* f8 - ff */
};

__constant__ unsigned int floodFillMasksSix[27]={
  0x20b,0x417,0x826,0x1059,0x20ba,0x4134,0x80c8,0x101d0,0x201a0,
  0x41601,0x82e02,0x104c04,0x20b208,0x417410,0x826820,0x1019040,0x203a080,0x4034100,
  0x2c0200,0x5c0400,0x980800,0x1641000,0x2e82000,0x4d04000,0x3208000,0x7410000,0x6820000
};

__constant__ unsigned int floodFillMasksTwentySix[27]={
  0x361b,0x7e3f,0x6c36,0x1b6db,0x3ffff,0x36db6,0x1b0d8,0x3f1f8,0x361b0,
  0x6c361b,0xfc7e3f,0xd86c36,0x36db6db,0x7ffffff,0x6db6db6,0x361b0d8,0x7e3f1f8,0x6c361b0,
  0x6c3600,0xfc7e00,0xd86c00,0x36db600,0x7fffe00,0x6db6c00,0x361b000,0x7e3f000,0x6c36000
};



__shared__ unsigned char sharedVoxelsX[2560];

/* utility kernels */

/* copy a strip of voxels, and the surrounding strips (if any), from device memory to shared memory */
__device__ void setupSharedVoxels(
    unsigned char *voxels,    /* input voxels in device memory */
    const int y,              /* y coordinate of the strip to be copied to shared memory */
    const int z,              /* z coordinate of the strip to be copied to shared memory */
    const unsigned int volDim /* lenght of the strip to be copied to shared memory */
    ) {

  /* Fill sharedVoxels with zeros */
  for (int i=0;i<10;i++){
    sharedVoxelsX[((volDim)*10*threadIdx.y)+10*threadIdx.x+i] = 0;
  }
  __syncthreads();

  /* Copy voxel values from voxels array to sharedVoxels */
  for (int i=0;i<3;i++) {
    for (int j=0;j<3;j++) {
      if (y+(j-1)>=0 && z+(i-1)>=0 && y+(j-1)<=(volDim-1) && z+(i-1)<=(volDim-1)) {
          sharedVoxelsX[((volDim+2)*9*threadIdx.y)+(i*3+j)*(volDim+2)+threadIdx.x+1]=voxels[(z+(i-1))*volDim*volDim+(y+(j-1))*volDim+threadIdx.x];
       }
    }
  }
}

__device__ unsigned int copySharedVoxelNeighbourhoodInt(const unsigned int volDim){
  unsigned int vnbh = 0;
  for (int i=0;i<3;i++) {
    for (int j=0;j<3;j++) {
      for (int k=0;k<3;k++) {
        if (sharedVoxelsX[(threadIdx.y*9*(volDim+2))+(j*3+k)*(volDim+2)+threadIdx.x+i]!=0) vnbh |= 1<<(i*9+j*3+k);
      }
    }
  }

  return vnbh; 
}

/* utility kernels */
__device__ unsigned int getViewAsIntZ(const unsigned int volDim,unsigned char label) {
  unsigned int view=0;
  int offset =(volDim+2)*9*threadIdx.y;
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,0);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,1);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,2);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x])) SET_BIT(view,3);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x])) SET_BIT(view,4);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,5);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,6);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,7);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,8);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,9);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,10);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,11);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x])) SET_BIT(view,12);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x])) SET_BIT(view,13);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,14);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,15);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,16);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,17);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,18);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,19);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,20);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x])) SET_BIT(view,21);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x])) SET_BIT(view,22);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,23);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,24);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,25);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,26);
  return view;
}

__device__ unsigned int getViewAsIntX(const unsigned int volDim,unsigned char label) {
  unsigned int view=0;
  int offset =(volDim+2)*9*threadIdx.y;
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,0);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,1);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,2);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,3);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,4);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,5);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,6);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,7);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x])) SET_BIT(view,8);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x])) SET_BIT(view,9);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x])) SET_BIT(view,10);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x])) SET_BIT(view,11);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x])) SET_BIT(view,12);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x])) SET_BIT(view,13);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,14);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,15);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,16);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,17);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,18);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,19);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,20);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,21);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,22);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,23);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,24);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,25);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,26);
  return view;
}

__device__ unsigned int getViewAsIntY(const unsigned int volDim,unsigned char label) {
  unsigned int view=0;
  int offset =(volDim+2)*9*threadIdx.y;
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,0);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,1);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,2);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+1)*(volDim+2)+threadIdx.x+0])) SET_BIT(view,3);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+0])) SET_BIT(view,4);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+0])) SET_BIT(view,5);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,6);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+1)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,7);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,8);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,9);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,10);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,11);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,12);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,13);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x])) SET_BIT(view,14);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,15);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+2)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,16);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,17);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,18);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,19);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,20);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(2*3+0)*(volDim+2)+threadIdx.x+0])) SET_BIT(view,21);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+0)*(volDim+2)+threadIdx.x+0])) SET_BIT(view,22);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+0])) SET_BIT(view,23);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,24);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(0*3+0)*(volDim+2)+threadIdx.x+2])) SET_BIT(view,25);
  if (HAS_LABEL(label,sharedVoxelsX[offset+(1*3+1)*(volDim+2)+threadIdx.x+1])) SET_BIT(view,26);
  return view;
}

__device__ int matchView2s_kernel(unsigned char *center,unsigned int view,unsigned int simpleview) {
  if (!GET_BIT(simpleview,8) || !GET_BIT(simpleview,26)) return 0;
  int t=(view&255)|((view>>9)&255);
  if ((topoTab[t][0]==1) && (topoTab[t][1]==1)) return 0;
  SET_2M_CRUCIAL(*center);
  if (topoTab[t][0]==0) {
    SET_SURF(*center);
  } else if (topoTab[t][1]>1) {
    SET_CURVE(*center);
  }
  return 1;
}

__device__ int match2s_kernel(const unsigned int volDim,unsigned char *center) {
  int ret=0;
  unsigned int view = getViewAsIntZ(volDim,S_ANY);
  unsigned int simpleView = getViewAsIntZ(volDim,S_SIMPLE);
  if (matchView2s_kernel(center,view,simpleView)) ret=1;
  view = getViewAsIntX(volDim,S_ANY);
  simpleView = getViewAsIntX(volDim,S_SIMPLE);
  if (matchView2s_kernel(center,view,simpleView)) ret=1;
  view = getViewAsIntY(volDim,S_ANY);
  simpleView = getViewAsIntY(volDim,S_SIMPLE);
  if (matchView2s_kernel(center,view,simpleView)) ret=1;
  return ret;
}

__device__ int matchView2_kernel(unsigned char *center, unsigned int view, unsigned int simpleView) {
  if (!GET_BIT(simpleView,8) || !GET_BIT(simpleView,26)) return 0;
  int t=(view&255)|((view>>9)&255);
  if ((topoTab[t][0]==1) && (topoTab[t][1]==1)) return 0;
  SET_2M_CRUCIAL(*center);
  return 1;
}

__device__ int match2_kernel(const unsigned int volDim,unsigned char *center) {
  int ret=0;
  unsigned int view = getViewAsIntZ(volDim,S_ANY);
  unsigned int simpleView = getViewAsIntZ(volDim,S_SIMPLE);
  if (matchView2_kernel(center,view,simpleView)) ret=1;
  view = getViewAsIntX(volDim,S_ANY);
  simpleView = getViewAsIntX(volDim,S_SIMPLE);
  if (matchView2_kernel(center,view,simpleView)) ret=1;
  view = getViewAsIntY(volDim,S_ANY);
  simpleView = getViewAsIntY(volDim,S_SIMPLE);
  if (matchView2_kernel(center,view,simpleView)) ret=1;
  return ret;
}

__device__ int matchView1_kernel(unsigned int view,unsigned int simpleview,unsigned int objectview,unsigned int crucialview) {
  unsigned int tmp1=((GET_BIT(view,2) && GET_BIT(view,4)) || (GET_BIT(view,3) && GET_BIT(view,26)));
  unsigned int tmp2=!(GET_BIT(objectview,2) && (!GET_BIT(simpleview,2) || GET_BIT(crucialview,2)));
  unsigned int tmp3=!(GET_BIT(objectview,3) && (!GET_BIT(simpleview,3) || GET_BIT(crucialview,3)));
  unsigned int tmp4=!(GET_BIT(objectview,4) && (!GET_BIT(simpleview,4) || GET_BIT(crucialview,4)));
  unsigned int tmp5=!(GET_BIT(objectview,26) && (!GET_BIT(simpleview,26) || GET_BIT(crucialview,26)));
  unsigned int tmp6=!(((view&0x723900)!=0) && (((view&0x3900)==0)||((view&0x720000)==0)));
  if (tmp1 && tmp2 && tmp3 && tmp4 && tmp5 && tmp6 && GET_BIT(view,26)) return 1;
  
  tmp1=((GET_BIT(view,2) && GET_BIT(view,0)) || (GET_BIT(view,1) && GET_BIT(view,26)));
  tmp2=!(GET_BIT(objectview,2) && (!GET_BIT(simpleview,2) || GET_BIT(crucialview,2)));
  tmp3=!(GET_BIT(objectview,1) && (!GET_BIT(simpleview,1) || GET_BIT(crucialview,1)));
  tmp4=!(GET_BIT(objectview,0) && (!GET_BIT(simpleview,0) || GET_BIT(crucialview,0)));
  tmp5=!(GET_BIT(objectview,26) && (!GET_BIT(simpleview,26) || GET_BIT(crucialview,26)));
  tmp6=!(((view&0x1e0f00)!=0) && (((view&0xf00)==0)||((view&0x1e0000)==0)));
  if (tmp1 && tmp2 && tmp3 && tmp4 && tmp5 && tmp6 && GET_BIT(view,26)) return 1;
  return 0;
}

__device__ int match1Z_kernel(const unsigned int volDim) {
  unsigned int view=getViewAsIntZ(volDim,S_ANY);
  unsigned int simpleView=getViewAsIntZ(volDim,S_SIMPLE);
  unsigned int crucialView=getViewAsIntZ(volDim,S_2M_CRUCIAL);
  unsigned int objectView=getViewAsIntZ(volDim,S_OBJECT);
  return matchView1_kernel(view,simpleView,objectView,crucialView);
}

__device__ int match1X_kernel(const unsigned int volDim) {
  unsigned int view=getViewAsIntX(volDim,S_ANY);
  unsigned int simpleView=getViewAsIntX(volDim,S_SIMPLE);
  unsigned int crucialView=getViewAsIntX(volDim,S_2M_CRUCIAL);
  unsigned int objectView=getViewAsIntX(volDim,S_OBJECT);
  return matchView1_kernel(view,simpleView,objectView,crucialView);
}

__device__ int match1Y_kernel(const unsigned int volDim) {
  unsigned int view=getViewAsIntY(volDim,S_ANY);
  unsigned int simpleView=getViewAsIntY(volDim,S_SIMPLE);
  unsigned int crucialView=getViewAsIntY(volDim,S_2M_CRUCIAL);
  unsigned int objectView=getViewAsIntY(volDim,S_OBJECT);
  return matchView1_kernel(view,simpleView,objectView,crucialView);
}

__device__ int match1_kernel(const unsigned int volDim) {
  unsigned int view = getViewAsIntZ(volDim,S_ANY);
  unsigned int simpleView = getViewAsIntZ(volDim,S_SIMPLE);
  unsigned int crucialView = getViewAsIntZ(volDim,S_2M_CRUCIAL);
  unsigned int objectView = getViewAsIntZ(volDim,S_OBJECT);
  if (matchView1_kernel(view,simpleView,objectView,crucialView)){
    return 1;
  }
  view = getViewAsIntX(volDim,S_ANY);
  simpleView = getViewAsIntX(volDim,S_SIMPLE);
  crucialView = getViewAsIntX(volDim,S_2M_CRUCIAL);
  objectView = getViewAsIntX(volDim,S_OBJECT);
  if (matchView1_kernel(view,simpleView,objectView,crucialView)){
    return 1;
  }
  view = getViewAsIntY(volDim,S_ANY);
  simpleView = getViewAsIntY(volDim,S_SIMPLE);
  crucialView = getViewAsIntY(volDim,S_2M_CRUCIAL);
  objectView = getViewAsIntY(volDim,S_OBJECT);
  if (matchView1_kernel(view,simpleView,objectView,crucialView)){
    return 1;
  }
  return 0;
}

__device__ int matchView0_kernel(unsigned char *center,unsigned int view,unsigned int simpleview,unsigned int
    crucialview2,unsigned int crucialview) {
  if (!((GET_BIT(view,12)&&GET_BIT(view,26)) || (GET_BIT(view,11)&&GET_BIT(view,4)) || (GET_BIT(view,13)&&GET_BIT(view,2)) || (GET_BIT(view,8)&&GET_BIT(view,3)) )) return 0;

  if (GET_BIT(view,12) && (!GET_BIT(simpleview,12) || GET_BIT(crucialview2,12) || GET_BIT(crucialview,12))) return 0;
  if (GET_BIT(view,26) && (!GET_BIT(simpleview,26) || GET_BIT(crucialview2,26) || GET_BIT(crucialview,26))) return 0;
  if (GET_BIT(view,11) && (!GET_BIT(simpleview,11) || GET_BIT(crucialview2,11) || GET_BIT(crucialview,11))) return 0;
  if (GET_BIT(view, 4) && (!GET_BIT(simpleview, 4) || GET_BIT(crucialview2, 4) || GET_BIT(crucialview, 4))) return 0;
  if (GET_BIT(view,13) && (!GET_BIT(simpleview,13) || GET_BIT(crucialview2,13) || GET_BIT(crucialview,13))) return 0;
  if (GET_BIT(view, 2) && (!GET_BIT(simpleview, 2) || GET_BIT(crucialview2, 2) || GET_BIT(crucialview, 2))) return 0;
  if (GET_BIT(view, 8) && (!GET_BIT(simpleview, 8) || GET_BIT(crucialview2, 8) || GET_BIT(crucialview, 8))) return 0;
  if (GET_BIT(view, 3) && (!GET_BIT(simpleview, 3) || GET_BIT(crucialview2, 3) || GET_BIT(crucialview, 3))) return 0;

  if (GET_BIT(view,26)) SET_0M_CRUCIAL(*center);
  return 1;
}

__device__ int match0_kernel(const unsigned int volDim,unsigned char *center){
  unsigned int view=getViewAsIntZ(volDim,S_ANY);
  unsigned int simpleView=getViewAsIntZ(volDim,S_SIMPLE);
  unsigned int crucialView2=getViewAsIntZ(volDim,S_2M_CRUCIAL);
  unsigned int crucialView=getViewAsIntZ(volDim,S_1M_CRUCIAL);
  if (matchView0_kernel(center,view,simpleView,crucialView2,crucialView)) return 1;
  return 0;
}

__device__ char getMostSignificantBit(unsigned int integer) {
  if (integer==0) return -1;
  unsigned char pos=0;
  if (integer>=1<<16) { integer>>=16; pos+=16; }
  if (integer>=1<< 8) { integer>>= 8; pos+= 8; }
  if (integer>=1<< 4) { integer>>= 4; pos+= 4; }
  if (integer>=1<< 2) { integer>>= 2; pos+= 2; }
  if (integer>=1<< 1) { pos+= 1; }
  return pos;
}

__device__ int hasOneConnectedTwentySixComponent(int vnbh){
  int out=0;
  if (vnbh==0) return out;
  unsigned int buffer=0;
  unsigned int filled=0;
  unsigned int visited=0;
  unsigned char center=GET_BIT(vnbh,13);
  UNSET_BIT(vnbh,13);
  unsigned char n=__ffs(vnbh)-1;
  SET_BIT(buffer,n);
  while (buffer>0) {
    n=__ffs(buffer)-1;
    SET_BIT(visited,n);
    filled|=vnbh&floodFillMasksTwentySix[n];
    buffer|=filled&~visited;
    UNSET_BIT(buffer,n);
  }
  if (filled==vnbh) out=1;
  else out=2;
  if (center==1) SET_BIT(vnbh,13);
  return out;
}

__device__ int hasOneConnectedSixComponent(int vnbh) {
  int out=0;
  if (vnbh==0) return out;
  unsigned int buffer=0;
  unsigned int filled=0;
  unsigned int visited=0;
  unsigned char center=GET_BIT(vnbh,13);
  UNSET_BIT(vnbh,13);
  unsigned char n=__ffs(vnbh)-1;
  SET_BIT(buffer,n);
  while (buffer>0) {
    n=__ffs(buffer)-1;
    SET_BIT(visited,n);
    filled|=vnbh&floodFillMasksSix[n];
    buffer|=filled&~visited;
    UNSET_BIT(buffer,n);
  }
  if (filled==vnbh) out=1;
  else out=2;
  if (center==1) SET_BIT(vnbh,13);
  return out;
}

/* 
 * find 18-neighbourhood, i.e. 6,2 geodesic neighbourhood
 */
__device__ int geodesic_neighbourhood(unsigned int nbhd){
  unsigned int tmp1 = nbhd&0x415410;
  unsigned int tmp2 = nbhd&0x2AA8AAA;
  unsigned int tmp3 = tmp1;
  if (tmp3>>4&1)  tmp3 |= 0xAA;
  if (tmp3>>10&1) tmp3 |= 0x80A02;
  if (tmp3>>12&1) tmp3 |= 0x208208;
  if (tmp3>>14&1) tmp3 |= 0x820820;
  if (tmp3>>16&1) tmp3 |= 0x2028080;
  if (tmp3>>22&1) tmp3 |= 0x2A80000;
  return tmp1|(tmp2&tmp3);
}

__global__ void findSimpleVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim, unsigned short *coordinatesList_dev) {
  int simpleTwentysix = 0;
  unsigned short y;
  unsigned short z;
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    y = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    z = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    setupSharedVoxels(voxels,y,z,volDim.x);
  }
    __syncthreads();

  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    unsigned int vneigh = 0;
    vneigh = copySharedVoxelNeighbourhoodInt(volDim.x);
    int connectedTwentysix=hasOneConnectedTwentySixComponent(vneigh);
    vneigh = ~vneigh;
    unsigned int geoNeigh = geodesic_neighbourhood(vneigh);
    int connectedSix=hasOneConnectedSixComponent(geoNeigh);
    simpleTwentysix=((connectedTwentysix==1) && (connectedSix==1));
 }
    __syncthreads();
    
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    if(IS_OBJECT(voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]) && simpleTwentysix){
        SET_SIMPLE(voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]);
    }
  }
}

__global__ void markSurfaceVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev){
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    unsigned short y = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned short z = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    setupSharedVoxels(voxels,y,z,volDim.x);

    __syncthreads();
    
    unsigned char center = sharedVoxelsX[((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1];

    if (IS_SIMPLE(sharedVoxelsX[((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1])) {
      if (match2s_kernel(volDim.x,&center)) {
        sharedVoxelsX[((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1] = center;
      }
    }
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[(threadIdx.y*9*(volDim.x+2))+4*(volDim.x+2)+threadIdx.x+1];
  }
}

__global__ void findSurfaceVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim, unsigned short *coordinatesList_dev) {
  unsigned short y;
  unsigned short z;
  unsigned int vneigh = 0; 
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    y = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    z = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    setupSharedVoxels(voxels,y,z,volDim.x);
  } 
    __syncthreads();
  
  if(blockIdx.x*2+threadIdx.y < numCoordinates){ 
    vneigh = copySharedVoxelNeighbourhoodInt(volDim.x);
    vneigh = ~vneigh;
  }  
    __syncthreads();

  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    if (IS_OBJECT(sharedVoxelsX[(threadIdx.y*9*(volDim.x+2))+4*(volDim.x+2)+threadIdx.x+1]) 
        && !(IS_SIMPLE(sharedVoxelsX[(threadIdx.y*9*(volDim.x+2))+4*(volDim.x+2)+threadIdx.x+1]))) {     
      unsigned int geoNeigh = geodesic_neighbourhood(vneigh);
      int connectedSix = hasOneConnectedSixComponent(geoNeigh);
      if (connectedSix>1) SET_SURF(sharedVoxelsX[((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1]);
      if (connectedSix==0) SET_SELECTED(sharedVoxelsX[((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1]);
    }
 }

    __syncthreads();

  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[(threadIdx.y*9*(volDim.x+2))+4*(volDim.x+2)+threadIdx.x+1];
  }
}

__global__ void unmarkNonSimpleVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim, unsigned short *coordinatesList_dev) {
  
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    unsigned short y = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned short z = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    setupSharedVoxels(voxels,y,z,volDim.x);

    __syncthreads();
    
    int numSelected=0;
    int numSurface=0;
    int offset = threadIdx.y*9*(volDim.x+2);
    UNSET_2M_CRUCIAL(sharedVoxelsX[offset+4*(volDim.x+2)+threadIdx.x+1]);

    for (int w=0;w<3;w++) {
      for (int v=0;v<3;v++) {
        for (int u=0;u<3;u++) {
          if (IS_SELECTED(sharedVoxelsX[offset+(w*3+v)*(volDim.x+2)+threadIdx.x+u])) numSelected++;
          if (IS_SURF(sharedVoxelsX[offset+(w*3+v)*(volDim.x+2)+threadIdx.x+u])) numSurface++;
        }
      }
    }
    if (numSelected==0 && numSurface==0) UNSET_SIMPLE(sharedVoxelsX[offset+4*(volDim.x+2)+threadIdx.x+1]);
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[(threadIdx.y*9*(volDim.x+2))+4*(volDim.x+2)+threadIdx.x+1];
  }
}

__global__ void find2CrucialVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev){
 
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    unsigned short y = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned short z = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    setupSharedVoxels(voxels,y,z,volDim.x);
    int offset = ((volDim.x+2)*9*threadIdx.y);
    const unsigned int sharedIndex = offset+4*(volDim.x+2)+threadIdx.x+1;
    __syncthreads();
    unsigned char center = sharedVoxelsX[sharedIndex];

    __syncthreads();

    if (IS_SIMPLE(sharedVoxelsX[sharedIndex])) {
      if (match2_kernel(volDim.x,&center)) {
        __syncthreads();
        sharedVoxelsX[sharedIndex] = center;
      }
    }
    __syncthreads(); 

    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[sharedIndex];
  }
}

__global__ void find1CrucialVoxelsX_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev) {
  if (blockIdx.x*2+threadIdx.y<numCoordinates) {
    unsigned int y=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned int z=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    unsigned int sharedIndex=((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1;
    setupSharedVoxels(voxels,y,z,volDim.x);
    __syncthreads();
    if (IS_SIMPLE(sharedVoxelsX[sharedIndex])) {
      if (match1X_kernel(volDim.x)) {
        SET_1M_CRUCIAL(sharedVoxelsX[sharedIndex]);
      }
    }
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[sharedIndex];
  }
}

__global__ void find1CrucialVoxelsY_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev) {
  if (blockIdx.x*2+threadIdx.y<numCoordinates) {
    unsigned int y=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned int z=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    unsigned int sharedIndex=((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1;
    setupSharedVoxels(voxels,y,z,volDim.x);
    __syncthreads();
    if (IS_SIMPLE(sharedVoxelsX[sharedIndex])) {
      if (match1Y_kernel(volDim.x)) {
        SET_1M_CRUCIAL(sharedVoxelsX[sharedIndex]);
      }
    }
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[sharedIndex];
  }
}

__global__ void find1CrucialVoxelsZ_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev) {
  if (blockIdx.x*2+threadIdx.y<numCoordinates) {
    unsigned int y=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned int z=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    unsigned int sharedIndex=((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1;
    setupSharedVoxels(voxels,y,z,volDim.x);
    __syncthreads();
    if (IS_SIMPLE(sharedVoxelsX[sharedIndex])) {
      if (match1Z_kernel(volDim.x)) {
        SET_1M_CRUCIAL(sharedVoxelsX[sharedIndex]);
      }
    }
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[sharedIndex];
  }
}

__global__ void find1CrucialVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev) {
  if (blockIdx.x*2+threadIdx.y<numCoordinates) {
    unsigned int y=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned int z=coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    unsigned int sharedIndex=((volDim.x+2)*9*threadIdx.y)+4*(volDim.x+2)+threadIdx.x+1;
    setupSharedVoxels(voxels,y,z,volDim.x);
    __syncthreads();
    if (IS_SIMPLE(sharedVoxelsX[sharedIndex])) {
      if (match1_kernel(volDim.x)) {
        SET_1M_CRUCIAL(sharedVoxelsX[sharedIndex]);
      }
    }
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[sharedIndex];
  }
}

__global__ void find0CrucialVoxels_kernel(const int numCoordinates,unsigned char *voxels,const dim3 volDim,unsigned short *coordinatesList_dev){
 
  if(blockIdx.x*2+threadIdx.y < numCoordinates){
    unsigned short y = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y];
    unsigned short z = coordinatesList_dev[4*blockIdx.x+2*threadIdx.y+1];
    setupSharedVoxels(voxels,y,z,volDim.x);
    
    __syncthreads();
    int offset = ((volDim.x+2)*9*threadIdx.y);
    unsigned char center = sharedVoxelsX[offset+4*(volDim.x+2)+threadIdx.x+1];
    __syncthreads();
    if (IS_SIMPLE(sharedVoxelsX[offset+4*(volDim.x+2)+threadIdx.x+1])) {
      if (match0_kernel(volDim.x,&center)) {
        sharedVoxelsX[offset+4*(volDim.x+2)+threadIdx.x+1] = center;
      }
    }
    __syncthreads(); 
    voxels[z*volDim.x*volDim.x+y*volDim.x+threadIdx.x]=sharedVoxelsX[offset+4*(volDim.x+2)+threadIdx.x+1];
  }
}

__global__ void findUnionOfSets(unsigned char *voxels,unsigned char *result,const dim3 volDim) {
  const unsigned int volIndex=blockIdx.x*blockDim.x+threadIdx.x;
  const unsigned int gatheredIndex=blockIdx.x;
  const unsigned int sharedIndex=threadIdx.x;
  __shared__ unsigned char sharedVoxels[256];
  __shared__ unsigned char sharedResult[256];
  sharedVoxels[sharedIndex]=voxels[volIndex];
  if ((sharedVoxels[sharedIndex] && !IS_SIMPLE(sharedVoxels[sharedIndex])) || IS_2M_CRUCIAL(sharedVoxels[sharedIndex]) ||
      IS_1M_CRUCIAL(sharedVoxels[sharedIndex]) || IS_0M_CRUCIAL(sharedVoxels[sharedIndex])) {
    sharedResult[sharedIndex]=1;
  } else {
    sharedResult[sharedIndex]=0;
  }
  if (sharedVoxels[sharedIndex] && !sharedResult[sharedIndex]) {
    sharedVoxels[sharedIndex]=0;
    sharedResult[sharedIndex]=1;
  } else {
    sharedResult[sharedIndex]=0;
  }
  if (sharedVoxels[sharedIndex]) sharedVoxels[sharedIndex]=S_OBJECT;
  __syncthreads();
  voxels[volIndex]=sharedVoxels[sharedIndex];
  //parallel reduction of result
  if (sharedIndex<128) sharedResult[sharedIndex]|=sharedResult[sharedIndex+128]; 
  __syncthreads();
  if (sharedIndex<64) sharedResult[sharedIndex]|=sharedResult[sharedIndex+64];
  __syncthreads();

  if (sharedIndex<32) {
    volatile unsigned char *smem=sharedResult;
    smem[sharedIndex]|=smem[sharedIndex+32];
    smem[sharedIndex]|=smem[sharedIndex+16];
    smem[sharedIndex]|=smem[sharedIndex+8];
    smem[sharedIndex]|=smem[sharedIndex+4];
    smem[sharedIndex]|=smem[sharedIndex+2];
    smem[sharedIndex]|=smem[sharedIndex+1];
  }
  if (sharedIndex==0) result[gatheredIndex]=sharedResult[0];
}

__global__ void gather_kernel(unsigned char *scattered,unsigned char *gathered) {
  const unsigned int volIndex=blockIdx.x*blockDim.x+threadIdx.x;
  const unsigned int gatheredIndex=blockIdx.x;
  const unsigned int sharedIndex=threadIdx.x;
  const unsigned int blockSize=blockDim.x;
  __shared__ unsigned int sharedVoxels1[256];
  sharedVoxels1[sharedIndex]=scattered[volIndex];
  __syncthreads();
  int stride=1;
  int multiple=2;
  while (multiple<=blockSize) {
    if (sharedIndex%multiple==0) {
      sharedVoxels1[sharedIndex]|=sharedVoxels1[sharedIndex+stride]; 
    }
    __syncthreads();
    stride*=2;
    multiple*=2;
  }
  __syncthreads();
  gathered[gatheredIndex]=sharedVoxels1[0];
}

__global__ void projectVolume_kernel(unsigned char *volume,dim3 volDim,unsigned char *yzPlane) {
  const unsigned int volIndex=blockIdx.x*blockDim.x+threadIdx.x;
  const unsigned int planeIndex=volIndex/volDim.x; 
  const unsigned int sharedIndex=threadIdx.x;
  __shared__ unsigned int sharedVoxels[256];
  sharedVoxels[sharedIndex]=volume[volIndex];
  __syncthreads();
  int stride=1;
  int multiple=2;
  while (multiple<=volDim.x) {
    if (sharedIndex%multiple==0) {
      sharedVoxels[sharedIndex]|=sharedVoxels[sharedIndex+stride];
    }
    __syncthreads();

    stride*=2;
    multiple*=2;
  }
  __syncthreads();
  if (sharedIndex%volDim.x==0) yzPlane[planeIndex]=sharedVoxels[sharedIndex];
}

__global__ void findStripCoordinates_kernel(unsigned char *yzPlane,dim3 planeDim,unsigned short *coordinateList,unsigned short *subListLengths) {
  const unsigned int index=blockIdx.x*blockDim.x+threadIdx.x;
  dim3 planeIndex;
  planeIndex.x=threadIdx.x%planeDim.x;
  planeIndex.y=(blockIdx.x*blockDim.x+threadIdx.x)/planeDim.x;
  const unsigned int listIndex=3*blockIdx.x*blockDim.x;
  const unsigned int sharedListIndex=threadIdx.x*3;
  const unsigned int subListIndex=blockIdx.x;
  __shared__ unsigned char sharedPlane[256];
  __shared__ unsigned int sharedCoordinateList[256*3];
  sharedPlane[threadIdx.x]=yzPlane[index];
  sharedCoordinateList[sharedListIndex]=0;
  if (sharedPlane[threadIdx.x]!=0) {
    sharedCoordinateList[sharedListIndex]=1;
    sharedCoordinateList[sharedListIndex+1]=planeIndex.x;
    sharedCoordinateList[sharedListIndex+2]=planeIndex.y;
  }
  __syncthreads();
  int stride=3;
  int multiple=2;
  while (multiple<=blockDim.x) {
    if (threadIdx.x%multiple==0) {
      int offset=2*sharedCoordinateList[sharedListIndex]+1;
      int count=sharedCoordinateList[sharedListIndex+stride];
      for (int i=0;i<count;i++) {
        sharedCoordinateList[sharedListIndex]++;
        sharedCoordinateList[sharedListIndex+offset+i*2]=sharedCoordinateList[sharedListIndex+stride+1+i*2];
        sharedCoordinateList[sharedListIndex+offset+i*2+1]=sharedCoordinateList[sharedListIndex+stride+i*2+2];
      }
    }
    stride*=2;
    multiple*=2;
    __syncthreads();
  }
  if (threadIdx.x==0) {
    subListLengths[subListIndex]=sharedCoordinateList[0];
  }
  __syncthreads();
  coordinateList[listIndex+threadIdx.x]=sharedCoordinateList[threadIdx.x];
  coordinateList[listIndex+blockDim.x+threadIdx.x]=sharedCoordinateList[blockDim.x+threadIdx.x];
  coordinateList[listIndex+blockDim.x+blockDim.x+threadIdx.x]=sharedCoordinateList[blockDim.x+blockDim.x+threadIdx.x];
}
